<<xAct`xTensor`

Define   manifold and metric

DefManifold[Σ,3,{i,j,k,l,α, β}];

DefMetric[1,S[-i,-j],DD,{"|","∂"}]

Define   constants and tensors

DefConstantSymbol/@{G,c,m,m1,m2,X1,X2,ν,r,rdot, γ,ω};

(*Relative separation vector*)
DefTensor[x[i],Σ]

(*Unit radial vector*)
DefTensor[n[i],Σ]

(*Tangential basis vector*)
DefTensor[λ[i],Σ];

(*Velocity and acceleration vectors*)
DefTensor[v[i],Σ]
DefTensor[a[i],Σ]

COM  definitions/ EOMs

massRules={m->m1+m2,X1->m1/m,X2->m2/m,ν->X1 X2};
gammaRule=γ->(G m)/(r (c^2));
revgammaRule=(G m)/(r (c^2))->γ;

(*Circular-orbit geometry*)
geometryRules={n[i] n[-i]->1,λ[i] λ[-i]->1,n[i] λ[-i]->0};
CircularOrbitRule={rdot->0};

(*Relative position*)
xRule=x[i]:>r n[i];

(*1PN orbital frequency*)(*Eq.(3.11)*)
omega1PNSqrule=(ω^2)-> (G m/(r^3))(1-(3-ν) γ);
omega1PNpow4rule=(ω^4)->((((G m/(r^3))(1-(3-ν) γ)))^2);
omega0PNSqrule=(ω^2)-> (G m/(r^3));
omega0PNpow4rule=(ω^4)-> (((G m/(r^3)))^2);

(*Velocity vector*)
velocityRule=v[i_]:>r ω λ[i];

(*Acceleration vector*)
accnRule=a[i_]:>-(ω^2) x[i];

vSq1PN=v[i] v[-i]/. velocityRule/.geometryRules/.omega1PNSqrule

v1PNSqRule=(v^2)->vSq1PN
v1PNpow4Rule=(v^4)->(vSq1PN^2)

Source  multipole  moments

DefTensor[STF2[i,j],Σ,Symmetric[{i,j}]];

DefTensor[STF3[i,j,k],Σ,Symmetric[{i,j,k}]];

(*Mass quadrupole (1PN)*)
Iij[-i,-j]:=m ν STF2[-i,-j] ( x[i] x[j] - (γ/42)(1+39ν) x[i] x[j] + (11/21)((r^2)/(c^2))v[i] v[j](1-3ν));

(*Mass octupole (0PN)*)
Iijk[-i,-j,-k]:=m ν(X2 - X1) STF3[-i,-j,-k](x[i] x[j] x[k]);

(*Current quadrupole (0PN)*)
(*Levi-Civita tensor*)
DefTensor[ϵ[-i,-j,-k],Σ,Antisymmetric[{-i,-j,-k}]];
Jij[-i,-j]:=m ν (X2 - X1)STF2[-i,-j]ϵ[-j,-α,-β](x[α]x[i]v[β]);

Derivatives

custom   derivative   rules

DefConstantSymbol[t];

td[a_+b_]:=td[a]+td[b];

td[a_ b_]:=td[a] b+a td[b];

td[a_^n_]:=n a^(n-1) td[a];

td[G]:=0;
td[m]:=0;
td[ν]:=0;
td[γ]:=0;
td[r]:=0;
td[c]:=0; 
td[X1] :=0;
td[X2] :=0;
td[x_?NumericQ]:=0;

(*STF projector treated as constant*)
td[STF2[__]]:=0;
td[STF3[___]]:=0;
td[ϵ[__]]:=0;

(*Orbital variables*)
td[x[i_]]:=v[i];
td[x[-i_]]:=v[-i];

td[v[i_]]:=a[i];
td[v[-i_]]:=a[-i];
td[a[i_]] := -(ω^2)v[i];

((Subscript[I, ij])^((2))) - 1PN

Iijdot=td[Iij[-i,-j]]

Iijdotdot = td[Iijdot] (* Perform substitutions using 1PN EOMs*)

Answer: m ν ({{STF2, ({{ ,  }, {i, j}})}}) (2 ({{v, ({{i}, { }})}}) ({{v, ({{j}, { }})}})+γ (-(23/21)+((9 ν)/7)) ({{v, ({{i}, { }})}}) ({{v, ({{j}, { }})}})-((2 G m ({{x, ({{i}, { }})}}) ({{x, ({{j}, { }})}}))/(r^3))+((G m γ ((149/21)-((23 ν)/7)) ({{x, ({{i}, { }})}}) ({{x, ({{j}, { }})}}))/(r^3)))

((Subscript[I, ij])^((3))) - 1PN

stfSymRule1=v[j] x[i]+v[i] x[j]:>2 x[i] v[j]

Iijdotdotdot = td[Iijdotdot] (* Perform substitutions using 1PN EOMs and truncate up to γ*)
(* Eq. C2a of BDI *)

((Subscript[I, ijk])^((1))) - 0PN

stfSymRule2=  x[i]x[j]v[k] + x[i]v[j] x[k] + v[i]x[j]x[k]:>3 x[i] x[j]v[k]

Iijkdot=td[Iijk[-i,-j,-k]]

((Subscript[I, ijk])^((2))) - 0PN

Iijkdotdot = td[Iijkdot]  (* We need up to 0PN only*)
(* Perform substitutions using 0PN EOMs *)

Answer : m (-X1+X2) ν ({{STF3, ({{ ,  ,  }, {i, j, k}})}}) (-((G m ({{x, ({{i}, { }})}}) ({{x, ({{j}, { }})}}) ({{x, ({{k}, { }})}}))/(r^3))+2 ({{v, ({{i}, { }})}}) (({{v, ({{k}, { }})}}) ({{x, ({{j}, { }})}})+({{v, ({{j}, { }})}}) ({{x, ({{k}, { }})}}))+({{x, ({{i}, { }})}}) (2 ({{v, ({{j}, { }})}}) ({{v, ({{k}, { }})}})-((2 G m ({{x, ({{j}, { }})}}) ({{x, ({{k}, { }})}}))/(r^3))))

((Subscript[I, ijk])^((3))) - 0PN

Iijkdotdotdot = td[Iijkdotdot] 
(* Perform substitutions using 0PN EOMs *)
(* Eq.  C1b of BDI*)

Answer: m (-X1+X2) ν ({{STF3, ({{ ,  ,  }, {i, j, k}})}}) (6 ({{v, ({{i}, { }})}}) ({{v, ({{j}, { }})}}) ({{v, ({{k}, { }})}})-((21 G m ({{v, ({{k}, { }})}}) ({{x, ({{i}, { }})}}) ({{x, ({{j}, { }})}}))/(r^3)))

((Subscript[I, ijk])^((4))) - 0PN

stfSymRule3=  x[i]v[j]v[k] + v[i]x[j] v[k] + v[i]v[j]x[k]:>3 v[i] v[j]x[k]

Iijkdotdotdotdot = td[Iijkdddot](* Perform substitutions using 0PN EOMs *)
(* Eq. C2b of BDI *)

Answer: m (-X1+X2) ν ({{STF3, ({{ ,  ,  }, {i, j, k}})}}) (-((60 G m ({{v, ({{i}, { }})}}) ({{v, ({{j}, { }})}}) ({{x, ({{k}, { }})}}))/(r^3))+((21 (G^2) (m^2) ({{x, ({{i}, { }})}}) ({{x, ({{j}, { }})}}) ({{x, ({{k}, { }})}}))/(r^6)))

((Subscript[J, ij])^((1)))

epsRule1=v[α] v[β] ϵ[-j,-α,-β]:>0
epsRule2=x[α] x[β] ϵ[-j,-α,-β]:>0

Jijdot = Collect[td[Jij[-i,-j]], ϵ[-j,-α,-β]]/. accnRule/. omega0PNSqrule

(* perform substitution and eps rules *)

Answer: m (-X1+X2) ν ({{STF2, ({{ ,  }, {i, j}})}}) ({{v, ({{i}, { }})}}) ({{v, ({{β}, { }})}}) ({{x, ({{α}, { }})}}) ({{ϵ, ({{ ,  ,  }, {j, α, β}})}})

((Subscript[J, ij])^((2))) - 0PN

Jij2 = td[Jij1](* Use 0PN EOM and eps rules*) (* C1f of BDI*)

((Subscript[J, ij])^((3))) - 0PN

Jijdotdotdot = td[Jij2](* Use 0PN EOM and eps rules*) (* Eq C1g of BDI *)

Answer: -((G (m^2) (-X1+X2) ν ({{STF2, ({{ ,  }, {i, j}})}}) ({{v, ({{i}, { }})}}) ({{v, ({{β}, { }})}}) ({{x, ({{α}, { }})}}) ({{ϵ, ({{ ,  ,  }, {j, α, β}})}}))/(r^3))

Energy Flux

(*Flux1PN = ((-G)/(c^5))((1/5)Iij3 * Iij3 + (1/(c^2))((1/189)Iijk4 *Iijk4 + (16/45)Jij3 * Jij3))*)

STFcontractRule1=STF2[-i,-j] STF2[-i,-j]x[i] v[j]x[i] v[j]:>(1/2) (r^2)(v^2)

Iij3Sq1PN = ( Iijdotdotdot * Iijdotdotdot/. STFcontractRule1//Expand)/.vSq1PNRule//Expand//Normal@Series[#,{γ,0,1}]&

((Answer: 32 (G^3) (m^5) (ν^2))/(r^5))+((64 (G^3) (m^5) γ (-106 (ν^2)+45 (ν^3)))/(21 (r^5)))

STFcontractRule3 = STF3[-i,-j,-k] STF3[-i,-j,-k] x[i] x[j] x[k] x[i] x[j] x[k]->(2/5) (r^6)
STFcontractRule4 = STF3[-i,-j,-k] STF3[-i,-j,-k] v[i] v[j] x[k] v[i] v[j] x[k]:>(2/5) (r^2) (v^4)
STFcontractRule5 = STF3[-i,-j,-k] STF3[-i,-j,-k] v[i] v[j] x[k] x[i] x[j] x[k]:>(2/5) (r^4) (v^2)

insideIijk4Sq = ((((Iijk4 * Iijk4/( (m^2) (((-X1+X2))^2) (ν^2)))//Expand)/. STFcontractRule3)/. STFcontractRule4)/. STFcontractRule5
(* substitute powers of v. Need only up to 0 PN so truncate powers of γ using Normal@Series *)

Iijk4Sq0PN = (((insideIijk4Sq/.vSq1PNRule/.vpow41PNRule//Expand)//Normal@Series[#,{γ,0,0}]&)//Expand)*( (m^2) (((-X1+X2))^2) (ν^2))

Answer: ((3042 (G^4) (m^6) (((-X1+X2))^2) (ν^2))/(5 (r^6)))